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ABSTRACT 


A time-variant  Wiener  filter  based  on  the  regionally  dispersive 
characteristics  of  long-period  surface  waves  is  presented.  The  dispersion 
curves  are  determined  with  maximum  entropy  spectral  analysis.  The  filter 
performance  is  tested  with  synthetic  chirp  waveforms  and  with  real  seismic 
data.  The  filter  enhances  the  estimation  of  signals  and,  in  particular,  the 
surface  wave  magnitude  measurability,  for  waveforms  with  signal-to-noise 
ratios  down  to  0 dB  RMS.  For  lower  signal-to-noise  ratios  the  signal  esti- 
mates are  less  reliable,  and  magnitude  estimates  may  be  biased.  Compared 
with  bandpass  filtering,  noise  rejection  improvement  ranges  from  3 to  9 db. 
The  filter's  ability  to  separate  multiple  signals  is  limited  by  filter  impulse 
response  side  lobe  interference.  The  filter  performance  is  furthermore 
limited  by  dispersion  of  noise,  by  non- stationary  noise,  by  the  filter's  inher- 
ent ability  to  generate  false  signals  from  broadband  noise  components,  by  the 
signal  bandwidth,  by  the  regional  dispersion  curve  variation,  by  narrowband 
filter  response  characteristics,  and  by  the  reliability  and  resolution  of  state- 
of-the-art  spectral  analysis  methods.  The  filter  appears  to  be  more  effective 
as  a signal  estimator  than  as  a detector. 
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INTRODUCTION 

Dispersion  of  surface  waves  has  been  used  in  signal  processing 
and  analysis  of  long-period  (LP)  seismic  waveforms.  For  instance,  correla- 
tion processes,  such  as  reference  waveform  and  chirp  waveform  matched 
filtering  (MF|,  compress  the  available  signal  energy  while  averaging  the  sup- 
:>•  sedly  random  noise,  to  yield  correlation  peaks  of  high  si gnal-to-noise  ratio 
S\h  i.  Because  of  signal  energy  compression  in  time,  these  techniques  do 
f r:  prove  the  estimate  of  the  original  signal  and  are  therefore  more  useful 
let*  i t on  purposes. 

This  report  presents  a method  to  enhance  the  estimate  of  an 
: - _ a . . t j y utilizing  prior  knowledge  of  its  dispersive  characteristics.  We 
. • s technique  dispersion-related  filtering  (DRF). 

The  most  elementary  form  of  DRF  consists  of  time-variant 
-.arrowr  and  filtering  along  the  signal's  presumed  dispersion  curve.  In  this 
• anner  considerably  more  noise  energy  can  be  rejected  than  in  stationary 
bandpass  filtering  of  the  waveform  over  the  entire  expected  signal  frequency 
band.  This  method  has  the  inherent  ability,  however,  to  generate  dispersed 
waveforms,  easily  mistaken  as  signals,  from  any  broadband  input,  including 
pure  noise.  For  this  reason,  DRF  is  not  well  suited  for  signal  detection;  its 
main  function  is  signal  estimation. 

For  a given  region- station  combination  the  expected  dispersion 
curve  may  be  obtained  by  overlaying  time-variant  signal  spectra  measured 
from  strong  events  occurring  in  a given  region.  The  spread  or  variance  of 
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the  dispersion  curve  ensemble  and  the  signal  bandwidth  at  each  point  in  time 
then  determine  the  bandwidth  which  should  be  applied  at  each  point  along  the 
dispersion  curve. 

The  signal  estimate  can  be  further  improved  by  performing 
time-variant  Wiener  filtering  (TVWF);  i.  e.  , for  each  point  along  the  disper- 
sion curve  the  expected  signal  power  and  the  expected  noise  power  are  bal- 
anced to  yield  a signal  estimate  of  minimum  mean  square  error.  In  particu- 
lar, TVWF  tends  to  reduce  signal  over- estimation  at  frequencies  of  relatively 
high  noise  power,  e.  g.  , at  0.  06  Hz  (17-second  micro-seismic  noise),  and  at 
frequencies  around  0.  02  Hz  (50  seconds). 

Besides  ambient  noise  rejection,  DRF  is  theoretically  capable 
of  separating  multiple  dispersed  signals,  provided  that  the  individual  disper- 
sion curves  can  be  resolved  by  spectral  analysis. 

It  is  clear  that  in  all  cases  the  DRF  and  TVWF  performance 
depends  largely  on  the  effectiveness  and  reliability  of  the  spectral  analysis 
method  applied.  Therefore,  part  of  this  study's  effort  was  dedicated  to  the 
use  of  a high-resolution  spectral  analysis  method,  the  maximum  entropy 
spectrum  (MES)  technique  (Burg,  1967). 

In  this  study  the  DRF  and  TVWF  methods  were  developed  and 
tested  on  synthetic  chirp  waveforms  and  on  beamed  waveforms  from  Sinkiang 
Province  seismic  events,  recorded  at  the  Alaskan  Long-Period  Array  (ALPA). 
The  use  of  this  data  base  permits  comparison  with  the  results  of  a previous 
matched  filtering  performance  study  (Unger,  1973).  The  emphasis  was 
placed  on  the  feasibility  of  dispersion-related  filtering  and  on  its  potential 
and  limitations  to  improve  the  estimation  of  long-period  seismic  signals  from 
noisy  waveforms. 

Discussion  of  the  above  topics  is  presented  in  the  following 
sections.  The  spectral  analysis  technique  and  its  results  are  discussed  in 
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Section  II.  Section  III  presents  the  DRF  and  TVWF  development  and  discusses 
various  design  aspects  and  performance  expectations.  The  DRF  and  TVWF 
performance  are  evaluated  in  Section  IV,  and  the  study  and  its  conclusions  are 
summarized  in  Section  V.  Related  literature  is  listed  in  Section  VI.  Some 
related  details  are  contained  in  the  appendices. 
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SECTION  II 
SPECTRAL  ANALYSIS 


A.  INTRODUCTION 

A prerequisite  for  the  success  of  dispersion-related  and 
Wiener  filtering  is  a reliable,  high-resolution  spectral  analysis.  The  maxi- 
mum entropy  spectrum  (MES),  first  presented  by  J.  P.  Burg  (Burg,  1967) 
and  comprehensively  discussed  by  T.  E.  Barnard  (Barnard,  1 975)  is  a high- 
resolution  spectrum  obtained  from  a relatively  small  number  of  lags  of  a 
waveform's  auto- correlation  function.  However,  the  reliability  and  the  reso- 
lution of  the  spectrum  depend  strongly  on  the  selected  input  parameters:  the 
waveform  gate  or  window  length  and  the  number  of  correlation  lags  used. 

For  time-variant  spectra  such  as  those  of  dispersed  wave- 
forms, the  choice  of  this  parameter  combination  is  even  more  delicate  than 
for  stationary  spectra.  On  the  one  hand,  a shortest  possible  gate  must  be 
selected  to  best  describe  the  spectral  variation  with  time.  On  the  other 
hand,  the  spectral  reliability  increases  with  the  ratio  of  gate  length  to  the 
number  of  lags  used,  while  the  resolution  improves  with  the  number  of  lags 
used. 

Yet  another  aspect  is  the  sampling  rate  of  the  waveform,  or 
equivalently,  its  Nyquist  frequency.  For  an  oversampled  waveform  the  MES 
algorithm  might  find  peaks  beyond  the  actual  signal  frequency  band.  This 
would  result  in  a loss  of  resolution  in  the  frequency  band  of  interest.  Thus, 
it  seems  appropriate  to  use  the  lowest  possible  sampling  rate  that,  accord- 
ing to  the  sampling  theorem,  still  adequately  describes  the  signal,  i.e., 
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1 / (2  W)  , 
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T = 
s 

where  T is  the  sampling  period  and  W is  the  signal  cut-off  frequency. 

ALPA  LP  signals,  for  instance,  usually  do  not  extend  (partly  due  to  the  instru- 
ment response)  beyond  0.  06  Hz,  so  that  we  could  low-pass  the  waveform  at 
0.06  Hz  and  sample  every  8 seconds  (0.  125  Hz  sampling  rate).  On  the  other 
hand,  over- sampling  provides  us  with  a greater  number  of  points  per  time 
unit,  which  might  permit  us  to  take  shorter  gates.  However,  the  fact  that 
these  points  are  no  longer  independent  may  nullify  this  assumed  advantage. 

Since  it  appears  to  be  difficult  to  analytically  establish  an  opti- 
mum MES  parameter  combination  (Barnard,  1975),  an  attempt  was  made  in 
this  study  to  empirically  arrive  at  some  best  set  of  parameters,  based  on 
spectral  analysis  of  synthetic  linear  chirp  waveform  combinations  and  real 
data.  The  "best"  parameter  combination  found  in  this  manner  was  then  ap- 
plied in  the  dispersion  curve  analysis  of  Sinkiang  Province  events.  The  noise 
spectral  analysis,  required  by  the  Wiener  filter,  demands  a different  set  of 
parameters. 

The  MES  parameter  analysis,  the  Sinkiang  region  dispersion 
curve  analysis,  and  the  noise  spectral  analysis,  respectively,  are  presented 
in  the  following  subsections.  The  maximum  entropy  spectra  used  in  the  anal- 
ysis were  computed  using  a variation  of  the  MES  algorithm,  the  Burg  tech- 
nique (Burg,  1968),  which  is  especially  efficient  for  the  computation  of  spectra 
from  short  waveform  windows.  It  bypasses  the  computation  of  the  actual 
auto-correlation  function  and  derives  the  spectrum  directly  from  the  coeffic- 
ients of  a waveform  prediction  error  filter  (Burg,  1 968;  Barnard,  1975).  The 
number  of  filter  coefficients  used  (the  filter  length)  equals  the  number  of 
auto-correlation  lags  in  the  original  MES  computation. 


II- 2 
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MAXIMUM  ENTROPY  SPECTRUM  PARAMETERS 


Combinations  of  linear  chirp  waveforms  were  synthesized  to 
test  the  reliability  and  the  resolution  of  the  maximum  entropy  spectra  (MES) 
and  to  obtain  a best  combination  of  window  length,  number  of  lags,  and  sample 
rate  as  parameters  for  the  MES  algorithm. 

Figures  II- 1 through  II-4  present  the  moving- window  maximum 
entropy  spectra  of  a 0. 01 5-0. 055  Hz  single  linear  chirp  waveform  and  of  a 
combination  of  two  waveforms,  for  various  parameter  combinations.  Each 
waveform  has  a 24-second  cosine  taper  at  both  ends.  For  each  point  in  time 
specified  (here,  every  8 seconds)  the  spectrum  is  computed  from  a waveform 
window  taken  symmetrically  about  that  point.  At  every  80  seconds  the  spec- 
trum is  fully  plotted,  with  its  power  density  in  the  direction  of  the  time  axis; 
the  power  density  reference  (origin)  is  at  the  corresponding  time  point  indi- 
cated,where  not  obscured  by  plotted  data,  by  a tickmark  on  the  upper  side  of 
the  time  axis.  The  spectra  are  held  constant  where  they  fall  below  their  refer- 
ence power  level  and  where  they  exceed  the  length  of  the  time  axis  to  the  far 
right.  The  reference  (origin)  power  level  (relative  to  1 computer  count)  anno- 
tated below  the  figure  depends  on  the  average  power  of  the  waveform;  the 
logarithmic  power  scale  (dB)  drawn  above  the  figure  depends  on  the  space 
between  the  plotted  spectra.  For  each  computed  spectrum  the  frequencies 
with  highest,  second  highest,  and  third  highest  power  density  are  indicated  by 
a circle,  triangle,  and  plus  sign,  respectively.  The  pattern  of  these  symbols 
establishes  the  basic  dispersion  curves;  the  individual  spectra  give  additional 
information  about  the  signal  bandwidth  along  the  dispersion  curves.  The 
actual  dispersion  curves  are  given  in  solid  lines. 

We  observe  from  Figure  II - 1 that  the  dispersion  of  a single 
linear  chirp  waveform  is  described  adequately  by  the  spectra  of  Figure  Il-lb; 
the  peaks  are  well  defined  and  narrow  and  accurately  follow  the  actual  disper- 
sion curve.  Increasing  the  number  of  lags  leads  to  spurious  peaks  as  show'll 
by  the  spectra  in  Figures  II-lc,  Il-ld,  and  Il-le.  In  Figure  II-lc  this  causes 
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MES  PARAMETER  ANALYSIS:  EFFECTS  OF  A HIGHER  SAMPLING  RATE 
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ambiguity  and  inaccuracy  in  the  description  of  the  dispersion  curve;  in  Fig- 
ures II- Id  and  II- le  it  suggests  the  presence  of  multiple  signals  with  parallel 
dispersions  in  addition  to  the  waveform  actually  present. 

Figure  II - 2 shows  the  spectra  for  a combination  of  two  wave- 
forms: the  same  linear  chirp  waveform  as  used  in  Figure  II- 1,  plus  a mono- 
chromatic waveform  starting  200  seconds  later.  We  notice  that  four  lags  are 
not  sufficient  to  describe  the  dual  dispersion.  The  spectra  obtained  with  six 
lags,  however,  adequately  describes  and  delineates  the  two  curves  and  espec- 
ially at  the  intersection  show  sufficient  resolution  power.  Again,  increasing 
the  number  of  lags  leads  to  spurious  spectral  peaks  beyond  those  that  cor- 
rectly follow  the  dispersion  curves,  with  Figure  II-2e  possibly  suggesting  the 
presence  of  a third,  weaker  waveform  dispersed  in  parallel  to  the  first  one. 

The  spectra  of  Figure  II- 3 were  obtained  using  shorter  wave- 
form window's  in  combination  with  the  six  lags  found  to  be  adequate  in  the 
previous  figure.  We  observe  that  shortening  the  gates  causes  more  scatter 
in  the  description  of  the  dispersion  curves. 

Next,  in  Figure  II-4,  we  investigate  the  effect  of  using  a higher 
sampling  rate.  Probably  because  of  the  fact  that  the  data  points  are  not 
independent,  this  does  not  allow  us  to  use  shorter  time  gates,  despite  the 
fact  that  a greater  number  of  waveform  points  is  used. 

Thus,  summarizing  the  effects  of  various  parameter  combin- 
ations on  the  spectra  of  known  waveform  combinations,  it  appears  thst  apply- 
ing the  Nyquist  sampling  rate  (8-second  period  or  0.  125-Hz  sampling  rate 
for  a 0.  06-Hz  signal  cut-off  frequency)  and  using  a 160-second  time  gate,  in 
combination  with  six  auto-correlation  lags  in  the  MES  computations,  ade- 
quately describe  and  separate  the  dispersion  curves  of  two  simultaneously 
occurring  signals.  For  a single  signal  only  four  lags  should  be  used.  It  is 
also  possible  that  the  dispersion  rate  (i.  e.  , frequency  change / signal  length) 
affects  the  MES.  This  was  not  investigated  in  this  study. 
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Now  the  problem  arises  that  with  real  signals  we  do  not  know 
how  many  spectral  peaks  are  authentically  present  at  a given  time,  so  that 
we  do  not  know  how  many  lags  to  use  to  adequately  describe  the  various,  and 
variant  number  of,  spectral  peaks.  For  instance,  LP  waveforms,  in  general, 
consist  of  multipath  signals,  may  contain  higher  modes,  and  are  affected  by 
noise  which  in  itself  may  possess  several  (possibly  non- stationary ) spectral 
peaks.  There  does  not  seem  to  be  a ready  solution  to  this  problem,  and  the 
best  that  one  can  do,  therefore,  is  to  experiment  with  the  parameters  and  to 
accept  the  most  plausible  outcome.  This  means  that  the  specification  of  region- 
al dispersion  curves  is  highly  subject  to  an  analyst's  spectral  interpretation. 

As  an  example,  we  will  investigate  the  Love  wave  beam  spectra 
of  one  of  the  Sinkiang  region's  strong  events  for  various  numbers  of  lags,  as 
given  in  Figure  II - 5 . To  a large  extent,  it  is  difficult  to  decide  which  spectral 
peaks  are  spurious  and  which  ones  are  not.  But  there  is  sufficient  consistency 
in  the  patterns  of  Figures  II-5c,  II-5d,  and  II-5e  to  suggest  that  at  least  two 
main  dispersion  branches  are  present.  Since  the  SNR  is  high,  about  30  dB 
RMS’",  the  noise  is  probably  of  minor  influence.  For  the  purpose  of  this 
study, this  analyst  settled  for  the  spectra  of  Figure  II-5c,  thus  assuming  that 
basically  two  dispersion  curves  are  simultaneously  present  for  at  least  part 
of  the  time. 

Accordingly,  the  8-second  sampling  time,  a 160-second  gate, 
and  six  auto-correlation  lags  will  be  used  in  the  further  spectral  analysis  of 
signals.  The  use  of  this  parameter  combination,  however,  may  cause  spur- 
ious peaks  to  appear  where,  in  reality,  only  a single  spectral  peak  is  present; 
and  it  may  cause  spectral  inaccuracies  where,  in  fact,  there  are  more  than 
two  peaks. 
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Since  the  noise  spectra  may  authentically  contain  more  than 
two  peaks,  the  noise  spectral  analysis  may  require  a different  set  of  param- 
eters. This  is  investigated  in  the  noise  spectral  analysis  subsection. 

C.  SINKIANG  REGION  DISPERSION  ANALYSIS 

MES  analysis  using  the  parameters  established  in  the  previous 
subsection  was  performed  on  the  Love  wave  beams  of  some  of  the  stronger 
events  of  the  Sinkiang  data  base  presented  in  Table  1 1 — 1 . This  data  base  is 
part  of  the  one  used  in  a previous  matched  filter  study  (Unger,  1 973).  The 
results  of  that  study  indicate  a good  correlation  between  the  signals  of  these 
events.  The  distance  between  event  epicenters  is  less  than  300  km  (Figure 
II- 6);  their  great  circle  epicentral  azimuths  toward  ALPA  differ  by  less  than 
2°.  The  event  numbers  used  in  the  matched  filter  study  have  been  adopted  in 
Table  II-  1 and  in  Figure  II - 6 . 

Figure  II- 7 shows  the  dispersion,  both  as  moving- window  spectra 
plots  and  as  group  velocity  curves  derived  from  these  spectra,  for  the  event 
SIN/170/17AL,  which  was  used  as  a reference  event  in  the  matched  filter 
study.  Figure  II - 8 presents  the  signal  spectra  of  four  other  relatively  strong 
events;  their  reference  power  levels  are  indicated.  The  figures  are  aligned 
with  respect  to  their  dispersion  curves.  Figure  II - 9 is  an  overlay  of  the  dis- 
persion curves  of  these  five  events.  Since  for  each  event  the  group  velocity 
scale  and  the  travel  time  are  slightly  different,  the  group  velocity  and  travel 
time  scales  of  event  SIN/170/17AL  were  taken  as  references  in  this  picture 
and  in  the  discussion  to  follow. 

We  observe  that  for  each  event  individually  the  dispersion  is 
rather  complex.  However,  over  most  of  the  signal  duration  - i.e.  , between 
roughly  4.  1 and  2.  9 km/sec  group  velocity  - the  five  events  show  consider- 
able consistency.  In  the  lower  frequency  part,  from  0.  016  Hz  to  0.  032  Hz 
(between  4.  1 and  3.  3 km/sec  group  velocity)  the  regional  dispersion  seems 
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well  defined  and  relatively  narrow  (less  than  0.  005  Hz  variation).  This  part 
is  followed  by  a discontinuity  near  3.4  km/sec  where  a branch  with  almost 
zero  dispersion  about  0.  040  Hz  occurs.  This  branch  terminates  approxi- 
mately where  a third  branch  (perhaps  a continuation  of  the  first  one)  of  higher 
frequencies  seems  to  start  up  (near  3.  1 km/sec).  At  the  3.4  km/sec  discon- 
tinuity, both  the  signal  bandwidth  and  the  spectral  uncertainty  (reflected  by 
the  amount  of  scatter  in  Figure  II - 9)  appear  to  have  increased.  For  each  in- 
dividual event  there  seems  to  be  additional  structure  beyond  the  dispersion 
curves  or  branches  described  here,  but  the  overlay  picture  does  not  show 
this  structure  to  be  very  consistent.  The  spectral  analysis  method  as  applied 
here  apparently  is  not  capable  of  resolving  possible  multiple  signals  in  the 
coda. 


An  extensive  discussion  on  the  possible  sei  smological  causes 
and  the  plausibility  of  the  measured  dispersion  structure  is  beyond  the  scope 
of  this  report.  For  this  study's  objective  - i.  e.  , to  demonstrate  the  feasi- 
bility of  dispersion- related  filtering  - the  main  feature  is  the  establishment 
of  a dispersion  pattern  which  to  a large  extent  is  regionally  consistent.  This 
same  feature,  of  course,  is  used  in  the  concept  of  matched  filtering.  In  this 
respect,  linear  chirp  waveform  matched  filtering  is  probably  less  successful 
than  reference  waveform  matched  filtering,  due  to  the  presence  of  more  than 
one  dispersion  branch. 

Since  the  "no-dispersion"  branch  about  0.  040  Hz  resembles 
part  of  an  oceanic  group  velocity  curve  (although  shifted  in  frequency.  Figure 
II- 10),  we  should  mention  the  fact  that  the  67°  great  circle  propagation  path 
includes  a 10°  to  15°  span  over  the  Arctic  Ocean. 
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FIGURE  II- 10 

LOVE  WAVE  DISPERSION  (OLIVER,  1962) 
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D. 


NOISE  SPECTRAL  ANALYSIS 


The  Wiener  filter  requires  the  expected  values  of  the  noise 
power  densities  present  at  each  point  during  signal  reception  (Section  III). 

One  way  of  obtaining  these  expected  values  is  to  assume  the  noise  to  be  sta- 
tionary over  an  interval  of  about  four  times  the  signal  duration,  and  to  meas- 
ure the  power  density  spectrum  of  a relatively  large  noise  gate  preceding  the 
signal.  Another  method  would  be  to  assume  the  noise  to  be  long  term  (at 
least  seasonally)  stationary,  and  to  average  the  spectra  over  an  ensemble  of 
typical  (seasonal)  noise  samples.  Because  of  time  limitations,  we  will  not  be 
able  to  verify  which  method  would  give  the  best  results,  and  we  will  proceed 
with  investigating  the  less  time  consuming  first  method. 

Thus,  we  will  examine  the  relatively  short  term  stationarity  of 
a typical  noise  sample.  Later,  in  Section  IV,  we  will  use  this  same  noise 
sample  in  the  actual  testing  of  the  dispersion- relation  filter.  First,  however, 
we  must  evaluate  which  parameters  to  use  in  the  MES  noise  analysis,  since  we 
do  not  know  a priori  how  many  distinct  spectral  peaks  will  authentically  be 
present  at  a given  time  in  the  noise  sample. 

Figure  II- 11  shows  the  spectral  analysis  for  various  MES  para- 
meter combinations,  of  the  first  2000  seconds  of  a typical  3800- second  ALPA 
horizontal  component,  north-looking  beam  noise  sample,  low-pass  filtered 
(with  a fourth-order  Butterworth  filter)  at  0.  06  Hz,  and  sampled  every  8 sec- 
onds. A plausible  result  seems  to  be  obtained  with  a 160- second  gate  and  a 
number  of  either  8,  10,  or  12  auto-correlation  lags.  These  patterns  show  a 
good  consistency;  yet  they  give  some  idea  of  the  degree  of  non-stationarity, 
under  the  (somewhat  ambiguous)  assumption  that  three  to  five  peaks  describe 
the  moving- window  spectra  more  or  less  adequately.  We  will  choose  ten 
lags  for  this  analysis. 
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FIGURE  II- 1 1 

NOISE  SPECTRA  FOR  VARIOUS  WES  PARAMETER  COMBINATIONS 

(PAGE  2 OF  2) 
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The  spectra  clearly  show  that  there  is  significant  variation, 
both  in  the  distribution  of  the  peaks  over  frequency,  as  well  as  in  the  peak 
magnitude  and  gate  average  power  levels.  Note  that  the  plotted  spectra  are 
independent;  i.  e.  , they  were  computed  from  non-overlapping  time  gates. 

Figure  11-12  is  the  MES  measured  from  the  first  1400  seconds 
(a  time  interval  corresponding  to  the  noise  gate  in  general  available  before 
the  first  P-wave  arrival  on  ALPA  records)  of  this  noise  sample.  This  spec- 
trum is  the  expected  power  density  spectrum  for  each  point  of  the  noise  occur- 
ring in  the  last  1800  seconds,  assuming  a signal  would  be  present  in  the  latter 
interval.  We  will  now  investigate  how  well  this  estimate  qualifies. 

Figure  11-13  shows  the  160-second  moving- window  noise  spectra 
of  the  second  part  of  the  noise  sample.  For  comparison,  the  expected  spec- 
trum has  been  drawn  in  as  a dotted  curve.  We  observe  power  density  devia- 
tions typically  on  the  order  of  6 dB  (occasionally  as  much  as  15  dB).  In  the 
Wiener  filter  design  discussion  in  Section  III,  we  will  see  how  these  deviations 
affect  the  filter  performance. 
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FIGURE  II-  I 3 

DEVIATION  OF  MOVING- WINDOW  SPECTRA  FROM 
EXPECTED  NOISE  SPECTRUM 
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SECTION  III 
FILTER  DESIGN 


A. 


INTRODUCTION 


This  section  presents  the  dispersion- related  filtering  (DRF) 
and  time-variant  Wiener  filtering  (TVWF)  developments.  Included  are  vari- 
ous design  aspects,  and  the  filter  performance  anticipation  based  on  the  re- 
sults of  the  spectral  analysis  obtained  in  the  previous  section.  Subsection  B 
discusses  the  basic  DRF  design;  Subsection  C presents  the  TVWF  process. 


B. 


DISPERSION- RELATED  FILTERING 


Dispersion-related  filtering  can  be  performed  as  a sampled, 
time-variant  convolution  along  the  signal's  presumed  dispersion  curve: 


N 


(t)  = ^ [x  (t)  * h (n,  t)]  * 6(t  - n • At) 

n=l 


(III- 1) 


where 


y(t)  is  the  DRF  output, 

x(t)  is  the  DRF  input, 

h(n,t)  is  the  DRF  impulse  response  at  time  n • At, 

§(t)  is  the  Dirac  function, 

N is  the  number  of  dispersion  curve  time  points, 

At  is  the  sampling  interval. 

In  our  DRF  implementation  the  convolution  is  performed  via  the  frequency 
domain: 
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F denotes  the  inverse  Fourier  transformation, 

X(jio)  is  the  Fourier  transform  of  X(t), 

H(n,oj)  is  the  (phaseless)  transfer  function  of  the  narrowband 
filter  applied  at  the  nth  dispersion  point. 


Thus,  the  DRF  consists  of  a bank  of  narrowband  filters  (NBFs), 
in  principle  one  NBF  for  each  dispersion  point,  and  the  DRF  output  is  synthe- 
sized from  adjacent  samples,  each  one  taken  from  the  output  of  its  corre- 
sponding NBF.  This  process  involves  one  forward  Fourier  transform  and  in 
principle  N inverse  transforms.  However,  at  points  where  the  NBF  speci- 
fication remains  unchanged  (for  instance,  if  part  of  the  dispersion  curve  is 
flat),  the  same  NBF  output  can  be  re-used  to  establish  the  next  DRF  output 
sample. 

The  shape  of  the  NBF  transfer  function  applied  is  basically 
rectangular;  its  gain  equals  one  for  the  frequencies  within  the  filter  band 
specified,  except  for  a 0.  002  Hz  cosine  taper  near  the  cut-off  frequencies. 

The  impulse  response  of  a rectangular  transfer  function  with  bandwidth  WHz 
has  a sin  x/x  envelope  (where  x = vW  t)  with  zeros  occurring  every  1 /W  sec- 
onds from  the  maximum.  The  cosine  taper  serves  to  reduce  the  amplitude  of 
the  impulse  response  sidelobe;  however,  it  also  changes  an  originally  speci- 
fied bandwidth  of  Hz  into  an  effective  bandwidth  of  approximately  W = 

- 0.002  Hz,  and  may  change  slightly  the  shape  of  the  impulse  response  also 
in  other  aspects. 

To  specify  the  desired  dispersion  band,  the  user  inputs  several 
dispersion  band  points  as  x-y  coordinates  in  inches,  based  on  his  measure- 
ment and  interpretation  of  the  moving- window  signal  spectra.  These  points  are 
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then  interpolated,  either  in  a linear  or  in  a quadratic  fashion  as  specified,  to 
establish  the  total  dispersion  band  to  be  used  by  the  DRF.  In  the  case  of  a 
linear  chirp  waveform,  a constant  bandwidth  of  linearly  increasing  center 
frequencies  can  be  applied.  However,  real  data,  at  least  as  shown  by  the 
Sinkiang  region  dispersion,  dictate  the  need  for  a dispersion  band  which  var- 
ies in  width  along  the  dispersion  curve. 

The  DRF  process  is  illustrated  with  the  filtering  of  a linear 
chirp  waveform  in  Figure  III  - 1 . Figure  III- la  presents  the  dispersion  band 
applied  in  the  DRF;  Figure  III  -lb  shows  the  original  waveform;  Figures  III-lc, 
III- Id,  and  III-le  are  the  respective  NBF  outputs  for  three  equally  spaced 
points  along  the  dispersion  curve;  Figure  III-  1 f is  the  DRF  output;  Figure 
III-  1 g is  the  error  trace  obtained  by  subtracting  the  signal  from  the  DRF  out- 
put. All  traces  are  plotted  on  the  same  scale.  Note  the  character  of  the  NBF 
outputs;  they  consist  of  a main  lobe  and  reduced  sidelobes;  the  node  interval 
(approximately  350  seconds)  corresponds  to  a 0.003  Hz  effective  (0,  005  Hz 
specified)  filter  bandwidth.  The  NBF  output  data  samples  which  synthesize 
the  DRF  output  are  indicated  with  arrows  and  dots.  A correction  factor, 
required  in  the  narrowband  filtering  of  time-variant  waveforms  (explained 
later  in  the  text  and  in  Appendix  A),  was  used  in  generating  the  NBF  outputs. 

We  observe  that  the  original  waveform  is  reproduced  fairly 
accurately;  there  is  some  distortion  due  to  NBF  end  effects.  The  first  and 
last  sample  of  the  DRF  output  waveform,  at  500  and  1500  seconds,  respec- 
tively, should  equal  zero,  but  the  NBF  is  not  infinitessimally  narrow,  thus 
permitting  the  signal  energy  at  neighboring  frequencies  to  pass  through  and 
to  establish  non-zero  values  at  those  points.  On  the  other  hand,  the  fact  that, 
for  the  NBFs  applied  over  approximately  the  first  and  last  100  seconds  of  the 
waveform,  the  signal  energy  of  the  neighboring  frequencies  on  one  side  is 
missing  results  in  signal  underestimation  over  those  intervals.  We  further- 
more notice  some  overshoot  (approximately  15%)  around  670  seconds  and 
1360  seconds,  and  a very  small  amount  of  ripple  over  the  rest  of  the  wave- 
form. 
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The  NBF  process  inherently  produces  amplitude  and  phase 
errors;  these  are  treated  in  more  detail  in  Appendix  A.  According  to  this 
treatment,  the  main  cause  of  the  NBF  amplitude  error  is  the  fact  that  the 
narrowband  frequency  energy  is  distributed  in  the  output  over  a relatively 
wide  NBF  response  interval.  This  error  can  be  corrected  by  applying  a 
bandwidth  and  dispersion  rate  dependent  correction  factor  (CORFAC).  Other 
effects,  which  could  not  be  compensated  for  at  this  point  in  the  development, 
cause  additional  amplitude  errors  as  well  as  phase  errors.  The  magnitude 
of  these  errors  seems  to  depend  on  the  nature  and  amplitude  of  the  input 
signal  and  on  the  bandwidth  of  the  NBF  applied.  In  Figure  III  - 1 the  amplitude 
error  (after  NBF  amplitude  correction)  is  less  than  5%,  excluding  the  over- 
shoot at  the  beginning  and  at  the  end  of  the  DRF  output.  The  relatively  large 
error  signal  displayed  in  Figure  III  - 1 g results  mainly  from  the  phase  error, 
leading  to  an  unrealistic  68.  6%  RMS  error  relative  to  the  RMS  signal.  Ap- 
pendix A shows  that  in  some  cases  the  uncor rectable  part  of  the  amplitude 
error  can  amount  to  30%,  i„e,,  Z.  3 dB  or  0.115  surface  wave  magnitude 
units.  This  potential  signal  distortion  must  be  weighed  against  the  signal 
estimate  improvement  obtained  from  the  filter's  noise  rejection  capability, 
which  ranges  approximately  from  3 to  9 dB  as  will  be  shown  next. 

At  each  point  along  the  dispersion  curve  the  amount  of  noise 
rejection  is  determined  by  the  width  of  the  dispersion  band  specified,  and  by 
the  actual  distribution  of  noise  energy  over  the  total  signal  frequency  band. 
For  instance,  the  Sinkiang  region  dispersion,  according  to  Figure  II - 9 , re- 
quires a bandwidth  of  at  least  0.  005  Hz  for  the  lower  frequency  part,  and  a 
bandwidth  of  possibly  as  much  as  0.  020  Hz  starting  at  the  dispersion  discon- 
tinuity at  3.4  km/sec  group  velocity.  With  respect  to  stationary  bandpass 
filtering  over  a 0.  015  to  0.  055  Hz  signal  band  (0.  040  H,.  Bandwidth),  and 
assuming  the  noise  to  be  white  over  that  frequency  band,  the  DRF  would  be 
capable  of  10  log  (0.040/0.005)  = 9 dB  noise  rejection  over  the  first  part  of 
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10  log  (0.040/0.020)  = 3 dB  over  the  second  part.  However,  some  signal 
energy,  contributed  by  frequencies  outside  the  dispersion  band  applied,  is 
possibly  sacrificed  in  the  filtering  process. 


These  amounts  of  DRF  noise  rejection  are  on  the  same  order 
as  that  of  a 9-element  beamsteer  filter,  and  as  the  matched  filter  SNR  gains 
over  bandpass  filtering  (Capon  et  al.  , 1969;  Strauss,  1973).  However,  it 
would  be  dangerous  to  use  the  DRF  as  a signal  detector  since  it  is  inherently 
capable  of  generating  false  alarms,  in  the  form  of  partial  chirp  waveforms, 
from  noise,  as  will  be  shown  in  the  filter  performance  analysis  (Section  IV). 
Also,  the  Sinkiang  signals  appear  to  have  their  highest  amplitudes  occur  dur- 
ing the  second  part  of  the  dispersion  curve,  where  the  DRF  noise  rejection  is 
less. 

To  apply  the  DRF  correctly  we  must  know  in  advance  the  signal 
start  time.  The  DRF  design  allows  for  "searching"  the  signal  by  sliding  the 
dispersion  band  through  the  received  waveform  data,  and  comparing  the  out- 
put peak  amplitudes  and  RMS  values  for  each  shift.  Supposedly,  the  shift 
yielding  the  highest  output  values  then  would  indicate  the  best  time  alignment 
between  the  signal  and  the  DRF  dispersion  band.  This,  however,  is  either  a 
computer  time,  or  a computer  core  consuming  process.  Moreover,  the  sig- 
nal spectrum  uncertainty  may  further  obscure  the  actual  signal  positions  in 
time.  The  procedure  works  w’ell  for  noise-free  linear  chirp  waveforms,  but 
loses  its  effect  under  marginal  noise  conditions  (about  0 dB  RMS  SNR).  The 
starting  problem  is  further  analyzed  in  Appendix  B where  alternative  attempts 
to  find  a signal's  start  time  also  are  described.  Among  these  alternatives 
are  the  time-domain  signal  phase  detection  method  and  the  MES  detection 
method.  All  these  methods  seem  to  have  their  limit  around  0 dB  RMS  SNR. 
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C.  TIME-VARIANT  WIENER  FILTERING 

The  DRF  reduces  the  error  in  the  estimation  of  signals  from 
LP  waveforms  by  eliminating  the  noise  energy  from  frequencies  outside  the 
signal  dispersion  band.  The  presence  of  noise  from  frequencies  inside  the 
dispersion  band,  however,  still  causes  some  interference,  in  general  leading 
to  signal  overestimation. 

A well-known  technique  of  reducing  the  estimation  error  is 
Wiener  filtering  (WF).  This  filtering  process  minimizes  the  mean-square 
estimation  error  by  balancing,  at  each  frequency,  the  signal  and  noise  power 
densities  (e.  g.  , Papoulis,  1967).  If  the  noise  is  uncorrelated  with  the  signal, 
the  time-variant  WF  transfer  function  is: 


H(t.w)  = 


Etf)g(t,  u) ) 


E <))q  (ti  to)  + E<f>N(t,w) 


(HI- 3) 


where 


H(t,o>)  is  the  time  variant,  phaseless  WF  transfer  function, 
E(i  (t,  a>)  is  the  expected  value  of  the  time-variant  signal  power 

O 

density  spectrum,  and 


(t,a>)  is  the  expected  value  of  the  time-variant  noise  power 
density  spectrum. 

We  can  also  write  the  WF  transfer  function  as: 


H(t,  ) = 


1 


1 + 


E<f>N(t,  w) 

(t,  oj) 


(HI- 4) 


To  show  the  effect  of  Wiener  filtering,  we  will  consider  some 
hypothetical  cases.  For  a noise-free  input  the  WF  transfer  function  equals 
one,  thus  passing  the  input  waveform  unchanged.  If,  at  a given  time,  the 
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expected  values  of  the  noise  and  signal  power  densities  are  equal  for  all  fre- 
quencies, the  transfer  function  equals  0,  5 for  all  frequencies;  i.  e.  , it  is  ex- 
pected that  in  the  mean,  for  equal  noise  and  signal  amplitudes,  the  best  esti- 
mate is  half  the  value  of  signal  plus  noise.  If  at  a given  time  and  for  a given 
frequency,  the  signal  power  density  equals  twice  the  noise  power  density,  then 
the  frequency  component  concerned  of  the  input  waveform  is  attenuated  with  a 
factor  0.  67.  In  general,  for  high  SNRs  the  filter  effect  is  little;  for  low  SNRs 
the  WF  performs  a frequency  and  time-dependent  attenuation. 

In  the  case  of  DRF  we  are  concerned  with  the  power  density 
ratios  within  a relatively  narrow  band  at  each  point  in  time  along  the  signal 
dispersion  curve.  For  a given  region  the  expected  values  of  the  time-varying 
power  densities  can  be  obtained  by  taking  the  moving- window  spectra  from  a 
strong  (reference)  event  signal,  or  by  averaging  the  spectra  of  strong,  re- 
gional event  signals.  The  expected  value  of  the  noise  power  densities  can  be 
obtained  by  assuming  the  noise  to  be  stationary  and  by  measuring  the  noise 
spectrum  of  a large  input  waveform  gate  just  prior  to  the  arrival  of  any  signal 
phase,  as  discussed  in  Section  II. 

Since  we  will  apply  this  filter  to  signals  of  different  strengths 
and  under  variant  noise  conditions,  we  must  adapt  the  spectral  ratio  term 
EVE*s  according  to  the  expected  SNR.  This  can  be  achieved  by  first  nor- 
malizing the  measured  noise  and  signal  spectra  to  unit  power,  i.  e.  , dividing 
the  power  densities  by  the  average  signal  and  noise  power,  respectively,  or 
by  the  squares  of  the  signal  and  noise  RMS  values.  We  then  multiply  this 
term  by  the  WF  design  parameter,  the  inverse  of  the  square  of  the  expected 
RMS  SNR: 


H (t,  cu ) = 


1 + 


ESNR 


PN(ClJ> 

Pg  ( t,oj ) 
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where 


P^((t>)  is  the  normalized,  supposedly  stationary  noise  power 
density  spectrum, 

P (t,w)  is  the  normalized,  time-variant  signal  power  density 
spectrum , 

ESNR  is  the  expected  RMS  SNR  of  the  input  waveform. 
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We  will  now  investigate  the  WF  effect  when  the  actual  power 
density  ratios  deviate  from  the  expected  ratios.  These  effects,  calculated 
from  Equation  III  - 5 , are  shown  in  Figure  III  - 2 for  various  true  SNR  situations. 
The  values  at  0 dB  error  are  the  correct  WF  attenuations  applied  to  the  input 
waveform  in  order  to  minimize  the  signal  estimation  mean-square  (m.  s,  ) 
error.  The  WF  gain  factor  becomes  more  erroneous  as  the  amount  of  power 
density  ratio  error  increases.  In  the  case  of  SNR  over e stimation,  the  WF 
still  reduces  further  the  DRF  signal  estimation  error.  For  SNR  underestima- 
tion, a point  may  be  reached  where  the  WF  attenuation  is  greater  than  twice 
the  correct  WF  attenuation,  indicated  by  the  vertical  dash-dot  line  at  0 dB 
SNR  error.  Beyond  that  point,  the  signal  estimation  error  caused  by  excess- 
ive WF  attenuation  is  greater  than  the  DRF  signal  overestimation  error.  The 
excessive  WF  attenuation  is  represented  by  the  dotted  parts  of  the  WF  gain 
curves.  Considering  this,  and  also  the  non-linearity  of  the  curves,  it  is  con- 
cluded that  less  error  is  made  when  overestimating  than  when  underestimating 
the  SNR. 

The  power  density  SNR  deviations  from  the  expected  SNR  at 
each  frequency  are  caused  by  two  factors: 

• The  moving-window  noise  spectrum  fluctuations  about  the  ex- 
pected value  of  the  noise  spectrum  (typically  6 dB,  occasionally 
10  to  15  dB). 

• The  mi sestimation  of  the  ESNR,  since  it  is  difficult  to  estimate 
the  SNR  from  low  SNR  waveforms  (probably  as  much  as  6 dB 
SNR  error  at  0 dB  true  SNR;  larger  errors  for  lower  true  SNRs). 
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Thus,  the  two  effects  combined  can  amount  to  deviations  spanning  the  range 
given  in  the  figure.  In  the  case  of  a 0 dB  true  SNR  at  a given  frequency,  this 
means  that  the  WF  gain  can  be  as  much  as  15  dB  (a  factor  5.  6)  too  low  or 
6 dB  (a  factor  2)  too  high  when  underestimating  and  overestimating,  respec- 
tively, the  power  density  SNR.  For  lower  true  SNR  situations,  these  errors 
are  larger.  For  a true  SNR  of  14  dB  or  higher  the  WF  has  little  effect,  and 
is  also  little  affected  by  SNR  estimation  errors. 

This  rapidly  puts  a limit  to  the  general  WF  performance;  at 
the  lower  SNRs,  where  improved  signal  estimation  is  desired  most,  the  ESNR 
is  difficult  to  estimate,  resulting  in  large  errors  in  the  WF  gain  factor  which 
in  turn  cause  larger  signal  estimation  errors.  From  the  curves  in  Figure 
III- 2 we  anticipate  unreliable  signal  estimation  for  waveforms  with  less  than 
0 dB  RMS  SNR. 

For  a given  waveform  the  SNR  is  probably  best  estimated  by 
visual  comparison  with  waveforms  obtained  through  simulation  of  various  SNR 
situations,  for  instance,  by  burying  the  regional  reference  signal  in  various 
levels  of  true  seismic  noise.  An  alternative  method  is  to  deduce  the  SNR 
from  the  bodywave  magnitude  and  measure  the  noise  level  just  prior  to  signal 
arrivals.  Because  of  and  m^-to-M^  conversion  uncertainties,  the  uncer- 
tainty in  the  ESNR  value  would  be  on  the  order  of  10  dB,  slightly  less  accurate 
than  SNR  estimation  directly  from  a waveform  of  0 dB  true  SNR,  but  probably 
equally  or  more  accurate  for  lower  true  SNR  waveforms. 

A third  method  is  to  compute  the  RMS  ratio  for  the  expected 
signal  and  noise  gates;  if  signal  and  noise  are  uncorrelated,  the  ESNR  can  be 
obtained  from: 


ESNR2  = RMSSIf 

rmsn 


(III-  6) 


III- 1 1 


where 


RMSSN  = signal-plus-noise  RMS  value,  and 
RMSN  = noise  RMS  value. 

Experiments  will  have  to  be  conducted  to  evaluate  this  method. 
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A.  INTRODUCTION 

In  this  section  we  will  evaluate  the  performance  of  both  the 
basic  DRF  concept  and  the  TVWF  process.  First,  the  DRF  and  TVWF  are 
applied  to  a linear  chirp  signal  with  added  seismic  noise  for  various  SNR 
levels.  Second,  we  will  apply  the  TVWF  to  a Sinkiang  region  reference  Love 
wave  beam  signal  under  various  noise  conditions,  and  to  the  waveforms  con- 
taining the  signals  of  other  Sinkiang  region  events.  In  the  third  subsection 
the  DRF  signal  separation  capability  is  tested  on  a combination  of  synthetic 
signals  and  tried  on  the  reference  signal. 

Since  the  Wiemr  filter  performance  is  greatly  affected  by  the 
power  SNR,  we  will  describe  the  test  conditions  in  terms  of  the  RMS  SNR, 
rather  than  the  peak-signal-to-RMS-noise  ratio  frequently  used  in  seismic 
analysis: 


RMS  SNR  = ZO  log,  „ gi*?  5'6  = 10  log,  „ 

10  RMS  noi  amp  10 


mean  sig  power 
mean  noi  power 


(IV- 1) 


For  seismic  signals,  roughly,  the  RMS  SNR  is  about  10  dB 
lower  than  the  peak-to-RMS  SNR  and  approximately  equals  the  peak-to-peak 
SNR.  For  uniform  envelope  single  chirp  signals  in  seismic  noise,  the  RMS 
SNR  is  3 dB  lower  than  the  peak-to-RMS  SNR  and  about  7 dB  higher  than  the 
peak-to-peak  SNR. 
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D. 


FILTER  PERFORMANCE  ON  LINEAR  CHIRP  WAVEFORMS 


Figure  IV- 1 shows  the  results  of  applying  the  basic  DRF 
(H  (o»,  t)  = 1,  followed  by  cosine  tapering,  within  each  NBF  passband)  to  wave- 
forms created  by  adding  seismic  noise  to  a synthetic  linear  chirp  waveform 
for  various  SNR  situations.  In  Figure  IV-2  the  DRF  outputs  from  the  low  SNR 
input  waveforms  are  plotted  on  a larger  scale  to  show  more  detail,  Down  to 
14  dB  RMS  SNR  the  DRF  rather  faithfully  reproduces  the  original  signal  with 
the  main  distortions  as  mentioned  in  Section  III:  waveform  tapering,  followed 
by  a 15%  overshoot  at  the  beginning,  and  the  same  phenomena  in  reverse  se- 
quence at  the  end  of  the  waveform;  some  ripple  along  the  rest  of  the  waveform; 
a relatively  larger  error  signal  due  to  phase  errors.  The  ripple  increases 
with  decreasing  SNR  due  to  a growing  interference  of  signal  and  noise.  From 
6 dB  SNR  down  the  distortion  increases  considerably,  in  ripple  as  well  as  in 
amplitude  error.  The  DRF  output,  however,  always  is  a much  better  approxi- 
mation of  the  original  signal  than  is  the  unfiltered  trace.  Notice  that  the  DRF 
produces  parts  of  a chirp  waveform  from  pure  noise;  if  there  is  no  information 
regarding  the  anticipated  signal  amplitude,  this  DRF  output  might  be  mistaken 
for  a signal. 

The  TVWF  performance  (using  the  transfer  function  given  by 
Equation  III- 5 , followed  by  cosine  tapering,  within  each  NBF  passband)  on  the 
same  input  waveform  is  displayed  in  Figures  IV- 3 and  IV-4.  In  each  SNR 
case,  the  true  RMS  SNR  was  input  as  the  ESNR  in  the  TVWF  algorithm  (see 
Section  III).  In  the  noise-only  case  a 0 dB  ESNR  was  used  to  test  the  filter 
output.  As  anticipated,  there  is  little  difference  in  performance  between  the 
DRF  and  the  TVWF  for  waveforms  of  14  dB  or  higher  SNR.  For  the  lower 
SNR  waveforms,  the  Wiener  filter  provides  a better  signal  estimate  than  does 
the  DRF;  over  the  entire  signal  duration  the  amplitude  error  is  smaller. 

Also  the  TVWF  apparently  can  generate  signal-like  waveforms  from  pure 
noise. 
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For  a quantitative  evaluation  of  filter  performance.  Figure 
IV- 5 displays  the  errors  in  maximum  amplitude,  excluding  the  overshoot  at 
the  beginning  and  end  of  the  signal,  for  the  unfiltered  trace,  the  0.  010-0,  060 
Hz  bandpass  filtered  (BPF)  trace,  the  DRF  output,  and  the  TVWF  output, 
respectively,  as  a function  of  RMS  SNR.  We  observe  that  the  TVWF  ampli- 
tude errors  are  in  good  agreement  with  the  errors  anticipated  based  on  the 
curves  of  Figure  1 1 1 - 2 and  on  the  amount  of  noise  spectral  variation  shown  in 
Figure  11-13.  This  amount  of  spectral  variation  establishes  the  limits  of  the 
TVWF  performance. 

The  graphs,  furthermore,  indicate  that  the  DRF  and  TVWF 
outputs  permit  accurate  magnitude  measurement  at  a significantly  lower  SNR 
than  is  possible  with  either  the  BPF  or  the  unfiltered  trace.  For  instance, 
measurement  of  the  peak  amplitude's  logarithm  (log  A)  within  _+ 0.  1 is  feasi- 
ble for  the  unfiltered  trace  only  for  input  waveforms  with  an  RMS  SNR  of  at 
least  16  dB:  with  TVWF  the  same  accuracy  of  measurement  can  be  obtained 
from  waveforms  with  an  RMS  SNR  as  low  as  5 dB.  A number  of  magnitude 
measurability  versus  SNR  relationships  are  listed  in  Table  IV-1.  This  table 
indicates  the  increase  in  magnitude  standard  deviation  caused  by  log  A mea- 
surement errors,  assuming  a distance  factor  uncertainty  (cr^)  of  0.3.  We 
notice  that  a log  A error  of  0.  1 does  not  significantly  increase  the  overall 
magnitude  error.  Magnitude  measurability  improvement  enhances  the  classi- 
fication of  seismic  events,  and  in  particular  should  improve  m^-M  discrimi- 
nation for  events  with  an  approximately  known  dispersion  curve. 

Comparing  the  merits  of  the  TVWF  with  alternative  methods, 
we  remark  that  possibly  similar  results  can  be  obtained  with  matched  filter- 
ing (MF),  since  the  amount  of  noise  rejection  is  on  the  same  order  as  that  of 
TVWF.  However,  the  large  variance  in  MF  gains  prevents  accurate  calibra- 
tion of  the  MF  output  (Strauss,  1 973),  and  the  MF  log  A measurement  accur- 
acy, therefore,  is  not  known.  The  MF  gain  probably  should  be  calibrated  by 
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Log  Peak  Amplitude  Error 


TABLE  IV-1 

SNR  REQUIRED  FOR  DESIRED  MAGNITUDE 
MEASUREMENT  ACCURACY*) 


M Accuracy  Log  A Minimum  RMS  SNR  Required  (dB) 

(assuming  cr^  = 0.  30)  Accuracy  Unfiltered  TVWF**) 


*)  For  linear  chirp  waveform  in  noise 
**)  SNR  presumed  known  exactly. 

***)  Extrapolated  value. 
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measuring  the  gain  when  a reference  waveform  is  ma#>  hed  to  itself,  rather 
than  the  present  procedure  of  calibration  with  regionally  average  MF  gains. 

The  above  error  analysis  is  based  on  the  combination  of  one 
simulated  chirp  waveform  and  only  one,  presumably  typical,  sample  of  seis- 
mic noise.  The  curves  in  Figure  IV-5,  therefore,  can  give  only  an  indication 
of  the  expected  filter  performance;  true  performance  characteristics  can  be 
obtained  only  from  analyzing  a number  of  noise  sample  and  signal  combina- 
tions, or,  perhaps  equivalently,  by  time- shifting  the  noise  waveform  with 
respect  to  the  signal.  However,  based  on  a quick  visual  inspection  of  other 
moving-window  noise  spectra,  we  expect  the  noise  spectral  power  density 
fluctuations  about  an  average  noise  spectrum  not  to  deviate  significantly  from 
those  displayed  in  Figure  11-13.  Since  the  peak  amplitude  errors  are  mainly 
determined  by  the  WF  mi sestimation  of  power  density  SNRs,  it  is  anticipated 
that  the  curves  in  Figure  IV-5  will  not  be  much  different  from  an  average 
TVWF  performance  curve  established  from  an  ensemble  of  signal  and  noise 
sample  combinations. 

Also,  this  TVWF  performance  is  based  on  an  exact,  a priori 
knowledge  of  the  ESNR,  which,  especially  for  low  SNR  waveforms,  is  a sen- 
sitive WF  design  parameter.  In  normal  filtering  situations,  additional  errors 
are  incurred  since  the  ESNR  must  be  estimated  by  visual  inspection  of  the 
input  trace,  or  by  deduction  from  the  bodywave  magnitude,  as  discussed  in 
Section  III.  Figures  IV-3  and  IV-6  indicate  that  the  SNR  of  a 6 dB  true  SNR 
waveform  can  probably  be  estimated  well  within  3 dB  accuracy.  Thus, 
according  to  Figure  III- 2 , the  WF  gain  error  due  to  SNR  misestimation  will 
be  small  (less  than  2 dB)  for  waveforms  of  6 dB  or  higher  true  RMS  SNR, 
but  they  may  become  significant  for  lower  SNR  waveforms. 

Besides  improving  the  surface  wave  magnitude  measurability, 
time-variant  Wiener  filtering  may  enhance  source  parameter  estimation. 

Love  wave  versus  Rayleigh  wave  energy  ratio  measurement,  and  possibly 
other  signal  analysis  and  classification  techniques. 
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TVWF  PERFORMANCE  ON  SEISMIC  WAVEFORMS 


Figure  IV- 6 shows  the  TVWF  performance  on  waveforms  com- 
posed of  the  Love  wave  beam  signal  of  the  Sinkiang  region  reference  event 
(Event  No.  3,  Table  II- 1,  Figure  II - 6 > and  added  seismic  noise  for  various 
SNR  situations.  Notice  the  shape  of  the  dispersion  band  applied,  made  to 
contour  the  regional  dispersion  variation  given  in  Figure  II - 9 ^ The  true  RMS 
SNR  value  was  input  as  the  ESNR  in  each  case. 
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We  observe  that  the  TVWF  reproduces  also  the  seismic  signals 
fairly  faithfully.  There  is  some  signal  distortion,  mainly  due  to  the  elimina- 
tion of  signal  energy  present  outside  the  filter  band  (see  Figure  II-8).  This 
again  shows  that  the  process  of  defining  the  signal  spectrum,  and  selecting 
the  corresponding  dispersion  band  to  be  used  by  the  filter,  is  somewhat  am- 
biguous; it  depends  on  the  analyst's  objective.  The  relatively  low  amplitude, 
and  for  the  most  part  non- sinusoidal,  error  signal  indicates  that  the  TVWF 
output  of  a natural  seismic  waveform  contains  little  phase  error.  For  the 
6 dB  and  lower  SNR  waveforms,  signal  overestimation  occurs  around  2000 
sec  travel  time.  Along  this  interval,  as  shown  by  the  corresponding  moving- 
window  noise  spectrum  at  2800  sec  in  Figure  11-13,  the  part  of  the  noise  spec- 
trum coinciding  with  that  of  the  signal  spectrum  (the  frequencies  around  0.  024 
Hz)  is  of  relatively  high  power,  and  is  about  6 dB  underestimated  by  the  ex- 
pected noise  spectrum.  According  to  Figure  III- 2,  this  results  in  signal 
overestimation  of  0.  5,  2,  and  6 dB,  respectively,  for  the  6,  0,  and  -6  dB 
SNR  waveforms,  in  good  agreement  with  the  amounts  that  can  be  measured 
from  the  TVWF  output  waveforms. 

The  errors  incurred  by  the  original  signal  peak  amplitude, 
i.e.,  the  amplitude  at  2300  sec  travel  time,  are  plotted  in  Figure  IV-7.  It 
appears  that  the  TVWF  performs  better  on  seismic  waveforms  than  on  the 
linear  chirp  waveforms.  This  is  in  part  due  to  a low  noise  level  relative  to 
other  parts  in  the  waveform  and  a rather  good  agreement  between  expected 
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and  actual  noise  spectrum  for  the  signal  peak  amplitude  time  and  frequency 
(3090  sec;  0.038  Hz,  Figure  11-13).  A different  time  alignment  between  the 
signal  and  the  noise  sample,  or  combinations  with  different  noise  samples, 
therefore,  could  result  in  somewhat  less  favorable  performance  curves. 
Nevertheless,  accurate  M^  measurement  (log  amplitude  measurement  within 
_+  0.  1 ) seems  possible  down  to  at  least  0 dB  RMS  SNR,  corresponding  approxi- 
mately to  surface  wave  magnitudes  reflecting  the  50%  detection  capability  of 
an  array  or  station. 

Using  the  same  regional  dispersion  band,  only  shifted  in  time 
to  align  with  the  various  signal  start  times,  the  TVWF  was  applied  to  the 
waveforms  containing  the  signals  of  the  other  events  listed  in  Table  II- 1.  The 
results  are  given  in  Figures  IV-8  and  IV-9.  For  each  event  the  expected 
noise  spectrum  and  its  average  power  level  were  measured  from  a 1440- sec 
gate  just  prior  to  the  first  P-wave  arrival.  The  expected  normalized  signal 
spectrum  is  determined  by  the  regional  dispersion  band,  by  the  power  densi- 
ties measured  and  interpolated  along  the  reference  signal  dispersion  curve, 
and  by  the  reference  signal  RMS  value,  as  described  in  Section  III.  The 
ESNRs  were  determined  by  visual  inspection  of  the  waveforms  and  compari- 
son with  the  simulated  SNR  conditions  of  Figure  IV-6.  The  signal  start  times 
were  estimated  with  the  techniques  described  in  Appendix  B. 

Since  we  do  not  know  the  actual  shape  of  the  signals,  but  antic- 
ipate them  to  be  similar  to  the  reference  signal  (Figure  IV-6),  we  can  merely 
observe  that  the  results  seem  very  plausible.  Certainly,  the  filtered  wave- 
forms give  the  impression  of  being  more  accessible  to  further  signal  analysis 
than  do  the  unfiltered  traces.  However,  we  do  not  know  if  and  how  much  sig- 
nal energy  from  frequencies  outside  the  dispersion  band  has  been  eliminated 
in  the  filtering  process.  We  notice  that  the  TVWF  outputs  for  the  Events  No. 
1,  2,  and  5,  which  are  within  20-25  km  from  the  reference  event  (Event  No. 

3,  Figure  II - 6 and  Table  II- 1)  are  indeed  very  similar  to  that  of  the  reference 
event.  The  events  at  100-300  km  distance  from  the  reference  event  display 
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slightly  different  signal  patterns.  Note  the  signal  similarity  between  Events 
No.  13  and  15.  The  signal  of  Event  No.  12  may  be  somewhat  uncertain  due 
to  the  low  SNR;  however,  in  shifting  the  TVWF  dispersion  band  in  time  through 
the  trace,  significant  signal  structure  in  the  TVWF  output  was  found  only 
within  the  expected  signal  time  gate  indicated  with  arrow's,  with  the  given  out- 
put yielding  the  highest  peak  amplitude.  Also,  the  signal  level  compares  well 
with  the  TVWF  output  of  the  -6  dB  SNR  waveform  in  Figure  IV-6.  The  filtered 
peak  amplitude  of  Event  No.  13  is  actually  larger  than  the  unfiltered  one,  pro- 
bably due  to  the  removal  of  interfering  signal  energy.  The  TVWF  output  of 
Event  No.  14  is  probably  distorted  due  to  unseparable  low  frequency  noise 
around  2100  sec  travel  time,  similar  to  the  situation  for  the  low  SNR  wave- 
forms in  Figure  IV-6. 

The  magnitude  changes  resulting  from  the  TVWF  are  summar- 
ized in  Figure  IV- 10.  No  conclusions  with  regard  to  magnitude  scatter  im- 
provement can  be  drawn  from  this  figure;  this  would  require  a larger  statisti- 
cal population  of  events. 

D.  SIGNAL  SEPARATION 

Figure  IV- 11  shows  the  DRF  separation  capability  for  a com- 
bination of  two  signals:  the  first  one  a linear  chirp  waveform,  the  second  one 
a monochromatic  signal.  The  DRF  adequately  resolves  the  two  signals  by 
separately  applying  the  two  corresponding  dispersion  bands  indicated  in  the 
figure.  The  chirp  signal  picks  up  part  of  the  monochromatic  signal  at  the 
intersection  of  the  dispersion  curves;  the  NBF  for  the  monochromatic  wave- 
form, however,  distributes  the  chirp  signal  energy  of  the  corresponding  fre- 
quency over  the  entire  time  interval  of  the  monochromatic  signal,  causing  a 
slight  increase  in  its  amplitude. 

Another  example  of  signal  separation  is  given  in  Figure  IV-12. 
Here,  the  input  consists  of  three  signals  with  parallel  linear  dispersion  curves, 
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simulating  multiple  signal  arrivals.  By  sliding  the  DRF  dispersion  band  for 
a single  signal  through  the  input  waveform,  each  of  the  three  signals  is  recov- 
ered whenever  the  DRF  dispersion  band  lines  up  with  a signal's  dispersion 
curve,  i.  e.  , at  300,  500,  and  700  sec,  respectively.  However,  there  is  a 
high  amount  of  distortion.  In  this  case,  the  separation  capability  is  severely 
limited  by  the  width  and  side  lobes  of  the  NBF  response  envelope  curve,  de- 
termined by  the  DRF  bandwidth.  On  the  one  hand,  we  want  to  make  the  width 
of  this  lobe  small  to  reduce  the  interference  caused  by  the  adjacent  main  lobes 
of  other  signals,  by  increasing  the  DRF  bandwidth.  On  the  other  hand,  this 
would  incur  both  stronger  side  lobes  which,  in  turn,  lead  to  interference,  and 
the  pick-up,  within  the  wider  DRF  band,  of  adjacent  signals  with  parallel  dis- 
persion curves.  Disregarding  the  effect  of  the  stronger  side  lobes,  we  can 
arrive  at  an  optimum  separation  bandwidth  Ws  ^ » yielding  a minimum  sepa- 
ration interval  Tg  . , as  follows.  The  separation  interval  Tg  is  greater 
than  the  lobe  width,  determined  by  the  DRF  bandwidth  W: 


T. 


2 W 


-1 


(IV- 2) 


But  also,  Tg  must  be  greater  than  an  interval  determined  by  W and  the  dis- 
persion rate  D,  as  indicated  in  the  figure: 
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The  above  is  only  true  if  the  DRF  band  is  sufficiently  narrow  to  preserve  the 
lobe  character  of  its  NBF  outputs.  According  to  Appendix  A,  this  condition 
seems  to  be  satisfied  for: 

W < D 1/2  . (IV -6) 

Thus,  the  optimum  bandwidth  cannot  satisfy  this  third  con- 
straint, leading  to  a compromise  between  losing  to  some  extent  the  single  lobe 

character  on  the  one  hand,  and  the  lobe  being  too  wide  on  the  other.  In  Figure 

- 5 

IV-12,  where  D = 4.  10  Hz/sec,  the  signals,  which  are  200  sec  apart,  were 
separated  with  an  effective  bandwidth  of  0.  008  Hz,  taking  into  account  the 
0.002  Hz  cosine  taper.  According  to  Equation  (IV-2),  this  results  in  a poten- 
tial separation  interval  of  250  sec,  which  still  satisfies  Equation  (IV-3): 

Tg  > 200  sec.  The  distortion  in  the  DRF  output,  therefore,  is  probably 
caused  mainly  by  the  interference  of  the  NBF  response  main  lobe  and  side 
lobes. 

Next,  Figure  IV- 13  shows  an  attempt  to  delineate  the  reference 
waveform  into  separate  signals,  more  or  less  similar  to  the  method  used  in 
Figure  IV- 11.  Moreover,  dispersion  band  1 is  applied  with  various  start 
times,  in  search  of  multiple  signals  with  parallel  dispersion  curves.  The 
chirp  signal  in  this  band  does  not  seem  to  extend  beyond  0.  035  Hz;  it  merely 
picks  up  signal  energy  from  band  2 at  the  intersection  of  the  two  dispersion 
curves,  and  some  other  scattered,  small  amounts  of  signal  energy.  Multiple 
signals  with  parallel  dispersion  curves  do  not  seem  to  be  present,  consistent 
with  the  spectral  analysis.  The  approximately  monochromatic  waveform  with 
a frequency  of  about  0.  040  Hz,  recovered  by  band  2,  is  also  consistent  with 
the  spectral  analysis.  Although  these  two  signals  seem  to  be  dominant  and 
together  probably  would  synthesize  most  of  the  input  waveform,  other  com- 
binations of  DRF  dispersion  bands  could  be  made  to  separate  or  recover  other 
signals  as  parts  of  the  composite  waveform. 
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A time- variant,  dispersion- related  Wiener  filter  (TVWF)  was 
developed  and  tested  on  synthetic  linear  chirp  waveforms  and  on  signals  from 
Sinkiang  Region  seismic  events.  The  maximum  entropy  spectrum  (MES)  tech- 
nique (Burg,  1968)  was  used  to  determine  the  Sinkiang  Region  dispersion 
curves. 


tions: 


These  experiments  led  to  the  following  conclusions  and  indica- 


•  The  MES  results  are  somewhat  ambiguous  and  depend  strongly 

on  the  MES  algorithm  parameters  used  (sample  rate,  wave- 
form gate  length,  number  of  prediction  error  filter  coeffic- 
ients). 


• The  definition  of  regional  dispersion  curves  is  subject  to  an 
analyst's  spectral  interpretation. 

• Despite  some  ambiguity,  the  MES  technique  provides  high- 
resolution  group  velocity  curves. 

• The  TVWF  is  effective  as  a signal  estimator  rather  than  as  a 
detector. 


• The  TVWF  enhances  the  estimation  of  signals  at  least  down  to 
0 dB  RMS  SNR.  In  particular,  it  considerably  improves  the 
measurability  of  surface  wave  magnitudes. 

• Below  0 dB  RMS  SNR,  the  estimates  may  become  unreliable 
due  to  noise  spectral  variation  with  time,  and  the  difficulty  of 
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estimating  the  waveform's  SNR,  which  is  a sensitive  parameter 
in  the  Wiener  filter  design. 

The  TVWF  noise  rejection  over  stationary  bandpass  filtering 
ranges  from  3 to  9 dB  depending  on  the  inherent  bandwidth  of  a 
signal  along  its  dispersion  curve. 

Dispersion- related  filtering  is  based  on  narrowband  filtering 
about  a known  dispersion  curve;  narrowband  filtering  of  time- 
variant  waveforms  produces  amplitude  and  phase  errors.  The 
amplitude  errors,  which  mainly  depend  on  the  filter  bandwidth 
and  the  dispersion  rate,  can  in  general  be  corrected  to  within 
1 dB;  in  unfavorable  cases  the  remaining  error  may  be  as  much 
as  2,  5 dB.  The  phase  error  could  not  be  corrected,  but  ap- 
pears to  be  small  for  natural  seismic  signals. 

The  filter's  separation  power  is  limited  by  the  widths  of  the 
filter  response  main  lobe  and  the  presence  of  side  lobes, 
determined  by  the  filter  bandwidth.  For  signals  with  parallel 
dispersion  curves,  the  minimum  separation  interval  and  the 
corresponding  bandwidth  are  determined  by  the  dispersion  rate. 
For  a 4*10  Hz/sec  dispersion  rate,  signals  separated  by  200 
sec  can  be  resolved  with  a bandwidth  of  0.  008  Hz;  however,  the 
output  contains  about  50%  amplitude  distortion  due  to  filter 
response  lobe  interference.  Signals  with  non-parallel  disper- 
sion curves  appear  to  be  better  separable. 

The  TVWF  requires  that  the  signal  start  time  be  known;  this 
may  be  found  with  any,  or  a combination,  of  the  following 
methods : 

deduction  from  given  source  location  and  time, 
sliding  the  TVWF  dispersion  band  through  the  waveform 
and  searching  for  the  maximum  RMS  output  value, 
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MES  analysis, 

instantaneous  signal  phase  detection. 

The  second  and  third  methods  appear  to  be  most  accurate,  but 
are  more  computer  time  and  core  consuming. 

• In  the  present  design  the  filtering  is  performed  in  the  frequency 
domain;  this  requires  in  principle  one  inverse  Fourier  trans- 
form for  every  dispersion  point  of  non-overlapping  bandwidths. 

For  a signal  with  a 1000-sec  dispersion,  sampled  at  2-sec 

-4 

intervals,  with  5*10  Hz  frequency  increments  and  a 0.  04  Hz 
bandwidth,  this  amounts  to  more  than  80  inverse  transforms. 
Alternatively,  one  may  conceive  the  filter  design  as  a time- 
domain  recursive  convolution,  using  a digital  resonant  filter 
technique;  this  method  should  be  considerably  faster. 

• The  TVWF  signal  enhancement  should  prove  useful  in  magnitude 
measurement,  Mg-m^  discrimination.  Love  wave  versus  Ray- 
leigh wave  energy  measurement,  source  parameter  studies, 
propagation  and  geological  structure  studies,  and  possibly 
other  signal  analysis  and  classification  techniques. 

• A statistical  filter  performance  evaluation  using  an  ensemble 
of  combinations  of  noise  samples  and  known  signals  is  required 
to  establish  the  full  range  of  filter  performance  characteristics. 
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B.  NARROWBAND  FILTER  AMPLITUDE  ERRORS 

Figure  A-l  shows  the  effects  of  narrowband  filtering  with  vari- 
ous bandwidths  along  the  dispersion  curve  of  a linear  chirp  waveform.  The 
narrowband  filter  (NBF)  transfer  function  and  its  response  characteristics 
were  described  in  Section  III.  The  shape  of  the  transfer  function  is  basically 
rectangular,  with  0.  002-Hz  cosine  tapers  near  the  cut-off  frequencies.  For 
sufficiently  narrow  filter  bands  (specified  later  in  the  text)  the  NBF  response 
to  a chirp  signal  input  has  a lobe  character  with  reduced  side  lobes;  the  node 
intervals  approximately  equal  the  inverse  of  the  effective  filter  bandwidth  (see 
also  Figure  III  - 1 ) . 

Figure  A- la  describes  the  variation  in  filter  bandwidth  applied 
along  the  dispersion  curve.  The  solid  lines  indicate  the  cut-off  frequencies, 
the  effective  filter  bandwidth  is  approximately  0.  002  Hz  smaller  due  to  the 
cosine  taper.  Figure  A-lb  is  the  input  signal;  Figures  A-lc  through  A-lh  are 
the  outputs  of  the  NBFs  applied  at  the  corresponding  points  along  the  disper- 
sion curve.  Figure  A-li  is  the  dispersion- related  filter  (DRF)  output 
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synthesized  from  samples  of  the  various  NBF  outputs.  Figure  A-lj  is  the 
error  signal  obtained  by  subtracting  the  input  signal  from  the  DRF  output. 

All  traces  are  plotted  on  the  same  scale. 

We  observe  that  the  wider  filter  bands  (0.  030-0.  020  Hz, applied 
from  500  to  850  seconds  in  the  time  trace)  reproduce  the  corresponding  signal 
peak  fairly  well.  The  convolution  of  the  chirp  signal  with  the  NBF  impulse 
response  creates  a small  amount  of  ripple  in  the  NBF  outputs.  From  0.  020 
through  0.  010  Hz  bandwidths  (applied  from  850  to  1300  seconds)  there  appears 
to  be  considerable  interference,  leading  to  amplitude  errors  of  up  to  30%,  as 
shown  in  the  DRF  output.  This  error  could  not  be  corrected. 

For  bandwidths  of  less  than  0.  009  Hz  with  the  given  dispersion 
rate,  the  NBF  outputs  assume  the  lobe  character  as  seen  in  Figure  III  - 1 . At 
this  point,  an  amplitude  error  is  created  due  to  the  fact  that  the  signal  energy 
carried  by  the  frequency  components  within  the  NBF  passband  is  distributed 
over  an  effective  output  time  interval,  determined  by  the  NBF  response  curve, 
which  is  greater  than  the  input  interval  established  by  the  time-frequency 
relationship  (the  dispersion  curve)  Over  the  passband  frequencies.  Based  on 
empirical  observations  and  on  theoretical  considerat  ons,  a formula  was  found 
which  satisfactorily  corrects  this  type  of  amplitude  error  for  monotonely  (not 
necessarily  linearly)  dispersed  waveforms.  To  correct  the  output  amplitude, 
it  is  multiplied  with  a correction  factor: 

COR  FAC  = W _1  D 1/2  , W < D 1/2 

= 1 , W > D 1/2  (A - 1 ) 


where 

D is  the  dispersion  rate  (Hz /sec)  over  the  NBF  passband  fre- 
quency components, 

W is  the  effective  NBF  bandwidth  (Hz). 
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This  formula  was  used  in  all  NBF  outputs  generated  in  this  study.  The  above 
conditions  are  indicated  in  Figure  A- la. 
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For  lobed  NBF  outputs,  this  correction  reduces  amplitude 
errors  to  less  than  5%,  not  counting  the  overshoot  at  the  end  of  a waveform 
as  reported  in  Section  III  and  also  presented  in  Figure  A-li„  The  formula 
was  tested  with  linear  chirp  waveforms  of  different  duration,  with  different 
dispersion  rates,  and  with  a cosine-modulated  linear  chirp  waveform. 

To  show  the  significance  of  this  correction,  we  calculate  the 

_5 

CORFAC  value  for  the  case  of  Figure  III  - 1 , where  D = 4*10  Hz /sec  and 
the  specified  bandwidth  of  0.  005  Hz  leads  to  an  effective  bandwidth  of  0.  003 
Hz.  The  result  is  a correction  factor  of  2.  1.  Thus,  without  the  correction, 
we  would  have  underestimated  the  signal  by  more  than  50%.  The  CORFAC 
values  are  plotted  as  a function  of  effective  filter  bandwidth,  for  several  dis- 
persion rates,  in  Figure  A-2. 

1 / 2 

The  condition  W < D has  a two-fold  significance.  First, 

from  Figure  A-l  it  is  observed  that  this  seems  to  be  the  condition  under 

which  the  NBF  output  has  the  lobed  character;  for  the  given  dispersion  rate 

(4’ 10  Hz/sec)  the  lobes  start  occurring  for  specified  bandwidths  of  less 

than  0.  009  Hz , cor  responding  to  an  effective  NBF  bandwidth  of  0.  007  Hz,  which 

1 /2 

is  only  slightly  higher  than  D = 0.0064  Hz.  For  wider  bands  the  lobed 
character  does  not  seem  to  be  present;  the  correction  factor  then  is  set  equal 
to  one  since  we  cannot  correct  for  the  amplitude  errors  occurring  in  that  case. 

1/2 

The  second  significance  of  the  condition  W < D is  that 

under  this  condition,  as  will  be  shown  shortly,  the  effective  output  duration 
is  greater  than  the  length  of  the  part  of  the  input  signal  comprised  by  the  NBF 
passband  frequency  components,  causing  the  output  amplitude  to  be  too  low 
with  respect  to  the  input  amplitude.  This  condition  requires  amplitude  cor- 
rection. 
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CORFAC  AS  A FUNCTION  OF  BANDWIDTH 
FOR  DIFFERENT  DISPERSION  RATES 


We  will  now  present  the  reasoning  which  led  to  Equation  (A-l). 
The  effective  input  duration  is  determined  by  the  NBF  bandwidth  and  the  dis- 


persion rate  over  this  band: 

T.  = W D . (A  - 2) 

We  define  the  effective  output  duration,  T , as  in  Figure  A-3.  It  is  the  inter- 

o 

val  over  which  a constant  amplitude  waveform  has  the  same  energy  as  the  lobed 
NBF  output  waveform;  the  constant  amplitude  equals  the  maximum  lobe  value. 

From  visual  inspection  of  the  lobed  waveforms  encountered,  it 
was  found  that  the  effective  duration  approximately  equals  one-half  of  the 
width  of  the  main  lobe,  which  is  2 W ^ . Therefore, 

T « W -1  (A-3) 

o 

was  adopted  as  a practical  value  for  the  effective  output  duration.  It  is  of 

interest  to  note  that  for  the  impulse  response  envelope  of  a rectangular  filter- 

band  (a  sin  x/x  function,  with  x = 7rW  t)  the  relation  is  exact:  T = W *. 

o 

For  any  other  lobed  curve,  Tq  can  be  calculated  from  numerical  integration. 
However,  using  the  above  approximate  value  returned  satisfactory  correc- 
tions, so  that  the  numerical  integration  was  not  performed  in  the  correction 
process.  It  now  follows  that: 

T > T.  for  W < D 1/2  . (A-4) 

o i 

The  correction  factor  then  follows  from  the  fact  that  the  input  signal  energy 

present  over  the  interval  T.  must  be  distributed  over  the  larger  interval  T , 

i o 

causing  the  output  amplitude  to  be  smaller  than  the  input  amplitude  by  a factor 

- 1 1/2 

(T.  Tq  ) . To  obtain  equal  input  and  output  amplitudes,  the  output  must 

be  multiplied  by: 


COR  FAC 


, W < D 


(A- 5) 
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T.  -V/2  = W-!D1/2 
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For  Tq  < T.  (i.  e.  , W>  D ) no  lohed  outputs  are  expected,  in  which  case 
we  cannot  perform  a correction  in  the  above  fashion,  and  the  correction  factor 
is  set  equal  to  one. 


In  the  present  stage  of  the  DRF  algorithm  Equation  (A-l)  was 
implemented  assuming  a linear  dispersion  curve  over  the  entire  signal  inter- 
val. The  dispersion  rate  is  automatically  determined  from  the  high  and  the 
low  frequencies  and  the  time  interval  of  the  dispersion  band  specified  and 
input  by  the  user.  This  method  can  easily  be  extended  for  the  case  of  non- 
linear dispersions  by  calculating  the  expected  instantaneous  dispersion  rate 
from  the  dispersion  band  center  frequencies.  All  data  presented  in  this  re- 
port were  processed  in  the  above  fashion;  for  the  measured  condition  W>D*^ 
the  correction  factor  is  automatically  set  equal  to  one. 


Thus,  the  correction  formula  was  effectively  applied  in  gener- 
ating the  DRF  and  time-variant  Wiener  filter  (TVWF)  outputs  of  a synthetic 
linear  chirp  waveform  (Figures  III-  1 , IV-1,  and  IV-3).  For  the  Sinkiang 
Region  seismic  waveforms,  the  dispersion  curve  variance  required  a band- 
width greater  than  the  square  root  of  the  average  dispersion  rate,  so  that  no 
correction  was  applied  (COR  FAC  = 1 ) in  that  case  (Figures  IV-6  and  IV-8). 


C.  NARROWBAND  FILTER  PHASE  ERRORS 

In  Section  HI,  Figure  III  - 1 g showed  a large  error  signal  despite 
good  amplitude  agreement  throughout  the  waveforms.  This  large  error  of 
rather  constant  amplitude  must  be  caused  by  a near- constant  phase  error,  as 
shown  below. 

Consider  a linearly  dispersed  waveform,  and  subtract  this  from 
the  same  waveform,  only  shifted  in  phase;  then  the  difference  signal  is: 


A-9 


* 


..  I 


M 


e(t)  = sin  (7rD  f + & ) - sin27rft 

= 2 sin  y-  cos  (ttD  f + y-  ) . (A-6) 

Thus,  a DRF  output  phase  shift  $ results  in  an  error  signal  of  constant  am- 
plitude , 

A =2  sin  ■—  , (A-7) 

e 

Tj-  (f) 

and  phase- shifted  an  amount  y + y with  respect  to  the  DRF  input. 

In  Figure  III  - 1 g the  error  signal  amplitude,  aside  from  end- 
effect  overshoot,  is  approximately  0.  7,  indicating  a phase  error  of  about  42  . 
An  overlay  of  the  DRF  output  trace  and  the  original  signal  trace  confirmed  a 
phase  error  of  this  order  of  magnitude. 

The  phase  difference  is  caused  by  the  convolution  of  the  chirp 
waveform  with  the  NBF  impulse  response.  This  convolution  is  rather  com- 
plicated; the  phase  difference  cannot  be  anticipated  and,  therefore,  cannot  be 
corrected  in  closed  form.  Notice  in  Figure  A-li  that  the  phase  error  seems 
to  be  small  for  the  relatively  wide  bands;  it  becomes  more  significant  for 
bandwidths  of  less  than  0.  020  Hz  (around  900  sec  on  the  time  axis). 

Probably  partly  due  to  the  wider  filter  bands  applied,  the  phase 
error  has  not  been  noticed  in  the  filter  outputs  of  actual  seismic  waveforms. 
Because  of  this  fact,  and  due  to  time  limitations,  the  phase  error  cause  and 
correction  were  not  pursued  further  in  this  study. 
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APPENDIX  B 


DETERMINING  THE  START  TIME 
OF  DISPERSED  SIGNALS  IN  NOISE 

A.  INTRODUCTION 

Operation  of  the  dispersion- related  filter  requires  prior  know- 
ledge of  the  start  of  a signal's  dispersion  curve.  Under  marginal  sign?.l-to- 
noise  ratio  (SNR)  conditions  this  may  be  difficult  to  establish.  Below,  four 
methods  of  determining  the  start  time  of  a signal's  dispersion  curve  under 
noisy  conditions  are  discussed. 

B.  DETERMINING  THE  START  TIME  FROM  TRAVEL  TIME  TABLES 

This  method  has  been  incorporated  in  the  EDIT  routine  of  the 
long-period  (LP)  signal  processing  package.  The  resulting  signal  start  times 
and  signal  duration  intervals  are  carried  in  the  LP  record  headers.  Signal 
start  and  stop  times  are  routinely  indicated  by  arrows  on  the  travel  time  axis. 
The  start  time  accuracy  is  on  the  order  of  100  sec  and  is  independent  of  the 
SNR. 

C.  DETERMINING  THE  START  TIME  WITH  THE  TVWF 

The  start  time  of  a signal's  a priori  known  or  assumed  disper- 
sion curve  can  be  found  by  sliding  the  TVWF  dispersion  band  through  the  wave- 
form. Correct  alignment  of  the  TVWF  dispersion  band  with  the  actual  signal 
dispersion  curve  should  yield  the  maximum  TVWF  output.  This  principle  is 
correct  for  signals  with  narrowly  defined  dispersion  curves.  In  the  case  of 
actual  seismic  signals  a relatively  wide  inherent  signal  bandwidth  and  the 
presence  of  unresolved  multiple  signals  may  cause  some  start  time  uncertainty. 


This  method  was  performed  on  both  a linear  chirp  waveform 

and  a Sinkiang  region  seismic  event  signal,  without  noise  and  with  a 0 dB 

RMS  SNR.  Figures  B-l  and  B-2  show  the  TVWF  RMS  output  as  a function  of 

presumed  signal  dispersion  start  time.  The  method  seems  to  work  well  even 

for  the  seismic  signal  in  a 0 dB  SNR  waveform;  the  start  time  for  the  chirp 

waveform  dispersion  is  estimated  correctly;  the  start  time  estimate  for  the 

seismic  signal  is  both  plausible  and  consistent.  The  plots  also  indicate  that 

a 50-second  timing  error  results  in  less  than  1 dB  amplitude  error,  i.e.  , 

less  than  0.  05  surface  wave  magnitude  error  (E  ). 

S 

However,  the  process  is  either  time  or  core  consuming  because 
of  the  large  number  of  inverse  Fourier  transforms  (one  for  each  different  NBF 
applied  along  the  dispersion  curve;  see  Section  III)  required  by  the  filter.  One 
has  the  choice  between  either  storing  the  relevant  intervals  of  all  NBF  outputs 
and  shifting  only  the  sampling  time,  or  repeating  the  entire  TVWF  process  for 
each  presumed  start  time.  The  first  method  requires  a high  amount  of  com- 
puter memory;  the  second  method  requires  considerable  processing  time. 

The  choice  depends  on  the  available  overall  processing  facilities  and  priorities. 

If  the  filter  were  to  be  designed  as  a recursive  digital  resonance 
filter,  the  processing  time  could  be  reduced  considerably,  due  to  the  possibil- 
ity of  fast  time-domain  convolution.  In  that  case  the  procedure  of  repeating 
the  TVWF  process  probably  is  preferable.  Moreover,  estimating  the  initial 
start  time  with  travel  time  tables  followed  by  a rough  probing  scheme  with  a 
subsequent  fine  search  could  be  scheduled  to  avoid  excessive  processing. 

The  results  in  Figures  B-l  and  B-2  are  only  an  indication  of 
the  start  time  accuracy  obtainable  with  this  method,  since  only  one  noise 
sample  was  used  in  conjunction  with  signals.  For  a statistical  evaluation  the 
method  should  be  tested  on  an  ensemble  of  noise  samples.  If  this  method 
shows  consistently  good  results,  also  for  lower  SNR  waveforms,  this  would 
enable  the  TVWF  to  be  used  also  as  a detector  rather  than  only  as  a signal 
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FIGURE  B-2 

DETERMINING  THE  DISPERSION  START  TIME  WITH  THE  TVWF; 
SIN/1 70/1 7AL.  (a)  SIGNAL  ONLY;  (b)  0 dB  RMS  SNR 
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estimator.  However,  as  it  is  possible  to  generate  false  signals  from  noise, 
it  may  also  be  possible  to  generate  false  start  times.  A statistical  evaluation 
will  have  to  indicate  the  potential  false  alarm  rate. 


D.  DETERMINING  THE  START  TIME  WITH  THE  MAXIMUM  ENTROPY 
SPECTRUM 

The  start  time  of  the  dispersion  curve  can  also  be  determined 
from  maximum  entropy  spectral  analysis  as  performed  in  Section  II.  Accord- 
ing to  Figure  B-3  the  high  resolution  of  this  technique  indicates  the  possibility 
of  picking  up  the  dispersion  curves  from  waveforms  of  0 dB  RMS  SNR,  prob- 
ably within  40  sec  accuracy.  Because  the  MES  has  to  be  performed  on  short 
overlapping  time  gates  which  slide  through  the  waveform,  also  this  method 
is  rather  time  consuming.  The  method  is  furthermore  sensitive  to  the  choice 
of  certain  MES  algorithm  parameters;  see  Section  II. 


E.  INSTANTANEOUS  ENVELOPE,  PHASE,  AND  FREQUENCY 
DETECTION 


A waveform  r(t)  can  be  expressed  as: 


r(t)  = R(t)  cos 
or,  equivalently,  as: 

r(t)  = R(t)  cos 


£ 2 7T  ^ f (t)  dt^ 

^ 2 7Tf0t  + j 


(B-l) 


(B-2) 


where,  see  Figure  B-4, 

R (t)  is  the  instantaneous  envelope  of  r(t), 

(|>(t)  is  the  instantaneous  phase  of  r(t)  with  respect  to  a 

monochromatic  waveform  of  frequency  fQ, 

f(t)  is  the  instantaneous  frequency  of  r(t), 

f is  an  arbitrary  reference  frequency, 
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FIGURE  B-3 


DISPERSION  START  TIME  DETERMINATION  WITH  THE  MES: 
(a)  LINEAR  CHIRP  WAVEFORM  IN  0 dB  RMS  SNR; 

(b)  SEISMIC  SIGNAL  IN  0 dB  RMS  SNR 
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The  envelope  and  phase  can  be  determined  with  the  waveform's 
Hilbert  transform  r(t),  a 90°  phase  - shift  operator  (e.  g.  , Papoulis,  1 965): 


r(t)  = R (t)  sin  [ 2 n fQt  + </>(t)J  . 


(B  - 3 ) 


by  the  relationships: 


R 


, [r2,t,  + 1/2 


'I' 


and 


r (t) 

<t>( t)  = arc  tan  — — - 27rf  t . 

r (t)  o 


(B-4) 


(B-5) 


The  frequency  is  then  found  by  differentiating  the  phase  function  after  modulo 
2tt  removal: 


«*>  - ^ ^ 


(B-6) 


For  a linearly  dispersed  waveform  the  frequency  is  a linear 

function  of  time.  Consequently,  the  phase  function  is  parabolic.  Since  the 

phase  cannot  exceed  an  interval  of  radians,  random  noise  with  its  frequency 

band  centered  at  f would  have  its  phase  fluctuations  confined  to  this  27r 
o 

radians  interval.  A systematic  phase  function  such  as  that  of  a dispersed 
waveform,  however,  can  be  "unwrapped"  by  eliminating  its  2n  modulo,  to 
expose  its  continuous  function.  This  may  result  in  a favorable  SNR  of  the 
phase  function,  depending  on  the  noise  spectrum  and  the  signal  dispersion 
curve.  The  phase  function,  therefore,  may  be  used  as  a detector,  and  could 
indicate  the  start  and  duration  of  LP  signals. 

In  reality,  however,  there  are  several  factors  that  may  reduce 
the  SNR  of  the  phase  function.  First,  the  center  frequencies  of  noise  and  sig- 
nal bands  may  not  coincide.  One  can  choose  f to  coincide  with  the  signal 
center  frequency;  in  that  case  the  offset  of  the  noise  center  frequency  causes 


the  noise  phase  function  to  have  a linear  bias.  Inversely,  one  may  choose 
to  be  the  noise  center  frequency;  then  the  signal  phase  will  have  a linear  bias. 
Experiments  will  have  to  indicate  which  choice  of  reference  frequency  is  best. 

Another  problem  is  that  the  noise  is  in  general  non-white,  due 
to  one  or  more  of  the  following  causes: 

• Propagated  noise  causing  dispersion. 

• Noise  coherence  due  to  instrument  response  and  beamsteering; 

in  particular,  possible  ringing  of  noise. 

• Natural  dominance  of  low-frequency  noise. 

This  causes  the  noise  phase  function  to  have,  over  short  intervals,  a some- 
what deterministic  character.  In  "unwrapping"  the  phase  the  peak  noise  phase 
values  then  may  become  higher  than  2t t,  thus  reducing  the  phase  SNR.  Pre- 
whitening will  not  give  a solution  since  this  would  affect  the  signal  spectrum. 

Figures  B-5  and  B-6  show  examples  of  envelope,  phase,  and 
frequency  detection.  After  making  the  phase  function  continuous,  some  of  the 
linear  bias  is  eliminated  by  taking  out  the  slope  between  the  first  and  the  last 
point  of  the  function.  The  frequency  is  derived  from  the  phase  function  before 
its  rotation.  Because  of  the  differentiation,  the  frequency  function  is  noisier 
than  the  phase  function.  In  particular,  it  shows  spikes  where  the  phase  func- 
tion has  discontinuities,  for  instance,  due  to  multiple  signals.  For  the  high 
SNR  waveforms,  the  parabolic  character  in  the  phase  function  can  be  clearly 
discerned,  despite  the  remaining  linear  bias  due  to  the  choice  of  f (0.040 
Hz);  the  frequency  function  is  clearly  linear.  Envelope,  phase,  and  frequency 
each  show  clearly  the  signal  onset  and  duration.  In  the  0 dB  SNR  waveforms, 
still  some  structure  can  be  recognized  in  each  of  these  functions,  but  although 
this  structure  is  probably  sufficient  to  call  a detection,  it  is  no  longer  possi- 
ble to  determine  the  signal  onset  and  duration.  Experiments  with  other  noise 
samples  yielded  somewhat  more  favorable  results;  the  merits  and  limitations 
of  this  method  will  have  to  be  determined  in  a separate  study. 
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SUMMARY 


Considering  the  above  analysis,  it  appears  that  determining  the 
dispersion  curve  start  time  is  probably  most  accurately  done  with  the  TVWF 
itself,  by  sliding  it  through  the  waveform  and  searching  for  the  maximum  out- 
put. This  method,  however,  is  time  or  core  consuming.  Next,  the  MES 
method  seems  to  give  the  best  results.  The  envelope,  phase,  frequency  de- 
tection method  needs  further  investigation.  The  travel  time  table  method  pro- 
vides initial  start  time  estimates.  The  signal  start  time  must  be  determined 
within  50  sec  accuracy  to  avoid  significant  magnitude  errors  due  to  dispersion 
curve  misalignment. 


